Beating patterns of filaments in viscoelastic fluids 
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Many swimming microorganisms, such as bacteria and sperm, use flexible flagella to move through 
viscoelastic media in their natural environments. In this paper we address the effects a viscoelastic 
fluid has on the motion and beating patterns of elastic filaments. We treat both a passive filament 
which is actuated at one end, and an active filament with bending forces arising from internal 
motors distributed along its length. We describe how viscoelasticity modifies the hydrodynamic 
forces exerted on the filaments, and how these modified forces affect the beating patterns. We show 
how high viscosity of purely viscous or viscoelastic solutions can lead to the experimentally observed 
beating patterns of sperm flagella, in which motion is concentrated at the distal end of the flagella. 



I. INTRODUCTION 

Eukaryotes use beating cilia and flagella to transport 
fluid and swim Cilia and flagella share a common 
structure, consisting of a core axoneme of nine doublet 
microtubules arranged around two inner microtubules. 
Molecular motors slide the microtubule doublets back 
and forth to generate the observed beating patterns. 
Cilia typically have an asymmetric beating pattern, with 
a power stoke in which the extended cilium pivots about 
its base, and a recovery stroke in which the cilium bends 
sharply as it returns to its position at the beginning of the 
cycle [2| . This stroke pattern is effective for moving fluid 
past the body of the cell, such as the epithelial cells of the 
airway or the surface of a swimming Paramecium. On the 
other hand, flagella typically have a symmetric planar or 
helical waveform [2| , such as the flagellum of bull sperm 
which exhibits a traveling wave with an amplitude that 
increases with distance from the head of the sperm ||. 
The shape of the centerline of a beating filament deter- 
mines the rate of transport or swimming. Many stud- 
ies show that this shape depends on the properties of 
the medium. For example, the flagella of human sperm 
have a slightly helical waveform in water. As viscosity 
increases via the addition of polymers, the waveform be- 
comes less helical and the amplitude flattens along most 
of the flagellum, with all the deflection taking place at 
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FIG. I: The waveform of a human sperm observed by Ishijima 
et al. [J in a solution with very high viscosity, 400 cP. 
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the free end [ij (Fig. [T]). Similar effects are observed in 
cervical mucus and other viscoelastic solutions [H, HJ ■ 

In this paper, we calculate the shape of a beating fila- 
ment as a function of the properties of the medium. Since 
the size scale for the filaments is tens of microns, and 
since typical swimming speeds are tens to hundreds of 
microns per second, the Reynolds number is very small, 
and inertia is unimportant. The medium is modeled 
with a single-relaxation time fading-memory model for 
a polymer solution 0]. We consider two different mod- 
els for the filament. First, we consider the passive case, 
in which one end of a flexible filament is waved up and 
down 0, H, We also consider a more realistic 

active model in which undulations arise from internal 
sliding forces distributed all along the length of the fil- 
ament [ill, Although active filaments are a closer 
representation of eukaryotic flagella, we show that many 
of the important features of swimming filaments in vis- 
coelastic media are also present in the simple passive 
filament. Furthermore, passive filaments may be more 
amenable to experiment and quantitative comparison 
with theory. Our main results are qualitative explana- 
tions for some of the beating shapes observed in sperm in 
high- viscosity and viscoelastic solutions [1, 0, H, EH, • 
These results indicate that the observed shape change 
with increasing viscosity is a physical rather than a be- 
havioral response. We reported some of our results for 
the case of an active filament in a viscoelastic medium 
with zero solvent viscosity in a previous publication [l5[ . 
In addition to presenting new results for the passive fil- 
ament in a viscoelastic medium, in the present paper we 
also systematically study the dependence of the shape of 
an active filament on solvent viscosity. 

To make the analysis as simple as possible, we work to 
linear order in the deflection of the filament away from a 
straight configuration. Although the power dissipated by 
a beating filament is second order in deflection, it is suf- 
ficient to calculate the shape to first order in deflection 
to accurately compute the power to second order. We 
calculate how this power varies as the beating pattern 
changes with increasing relaxation time. There are also 
important changes in the swimming velocity due to shape 
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change. Even in a linearly viscoelastic fluid, for which the 
swimming speed of a filament with prescribed beating 
pattern is the same as in a Newtonian fluid [IB] , changes 
to the beating patterns due to the viscoelastic response 
can lead to changes in swimming velocities. However, as 
has been previously emphasized [HI, [T3|, in viscoelas- 
tic fluids the swimming velocity also receives corrections 
from the nonlinearity of the fluid constitutive relation at 
second order, and both these nonlinear corrections and 
the changes to the beating patterns must be taken into 
account to correctly calculate the swimming velocity [l5[ . 

We begin our analysis by reviewing commonly used 
theories for the internal forces acting on passive and ac- 
tive filaments. Then we introduce the Oldroyd-B fluid, 
the fading-memory model we will use throughout this 
paper. After providing context with a brief discussion 
of resistive force theory for Newtonian viscous fluids, we 
discuss resistive force theory for small amplitude motions 
of a slender filament in an Oldroyd-B fluid. Using the 
balance of internal and external forces, we calculate the 
shape of a beating filament for both the passive and ac- 
tive cases. We close with a discussion of the limitations 
of our analysis and the implications of our results. 



II. MODELS 

A. Internal forces acting on flexible filaments 

The internal forces acting on a filament may be spec- 
ified either by the energy as a function of the curvature 
or a constitutive relation for the moment acting on a 
cross-section of the filament. The two approaches are 
equivalent, and in this paper we use the moment on a 
cross-section. 



1. Passive Filament 

We first consider the passive elastic filament. Let r(s) 
denote the path of the centerline of an inextensible fila- 
ment of length L. Then the moment M due to internal 
stresses acting on the cross-section of a bent filament is 



M = Ay' x r" 



(1) 



where A is the bending modulus and primes denote dif- 
ferentiation with respect to s [Lsl ]. Since we only con- 
sider planar filaments in this paper, it is convenient to in- 
troduce the right-handed orthonormal frame {ei,e2,z}, 
where ei = r' is the tangent vector of r(s), and z is 
the outward normal to the plane of the filament. Thus, 
K = 62 • is the curvature of the the curve r(s). Note 
that the sign of k is meaningful since &i does not flip 
at inflection points. Throughout this paper, we use a 
coordinate system for which the straight filament lies in 
the direction, and consider only planar motion of the 
filament in the e x -e y plane. Further, we restrict our- 
selves to small displacements of the filament, so that 
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FIG. 2: The dotted line is the centerline of the flagellum when 
it is straight; the circles divide it in quarters. The solid lines 
represent the bent centerline and r±; the dots divide each 
line in quarters as measured along the respective contours. 
All three curves have length L. 



r(s) « se x + h(s)e y . For such small displacements, we 
may approximate k « h" and M ks zAh". The elastic 
force per unit length on the filament is found by moment 
balance, 



M' +ei x F = 0, 



(2) 



where F is the force acting on a cross-section at s. For 
small displacements, the force per unit length due to in- 
ternal stresses is therefore fj nt = F' ~ —Ah""e y . Eq. @ 
also allows a tangential force that may be identified as a 
tension. For small displacements, the tension force is 
higher order in displacements than the normal forces, 
and so we ignore it in this work [Tgj | . In the small 
Reynolds number flows associated with swimming cells, 
the total force on each element of the filament must be 
zero, f cxt + fint = 0. As we shall see below, to deter- 
mine the shape of the filament, it is sufficient to consider 
external forces which are linear in the velocity of fila- 
ment elements. This conclusion is true even in nonlin- 
early viscoelastic fluids. Therefore force balance yields a 
fourth order linear differential equation to be solved for 
the time-dependent filament shape. 

The differential equation must be supplemented with 
appropriate boundary conditions. For the free end, we 
have 



Ah"'(L) = 
Ah"{L) = 0, 



(3) 
(4) 



corresponding to zero force and torque, respectively. For 
the driven end, we may consider prescribed angle, 



or prescribed force 



h'{0) = ecos(wi) 
h(Q) = 0, 

-Ah"'(0) = ecos(wf) 
h(0) = 0. 

2. Active Filament 



(5) 
(6) 

(7) 
(8) 



For a simple model of a sperm flagellum which incorpo- 
rates active bending forces, we consider a swimmer con- 
sisting of a flagellum which always lies in the plane, and 



3 



disregard the presence of a head. We model the flagellum 
as two inextensible filaments of length L which are paral- 
lel and have a fixed, uniform separation a. The filaments 
can slide past each other, which causes the flagellum to 
bend (Fig. [2]). Define the paths of the two filaments as 

r± =r±(a/2)e 2 , (9) 

where the frame {ei, 62, z} again applies to the centcrlinc 
r(s) of the filament (Fig. [2]). Equation ^ implies ds + — 
ds_ = — ands, where s±(s) is the arclength of r±(s). For 
example, n < in the interval 0<s<L/2of the curve r 
in Fig.O Thus, the arclength s + {L/2) of r + (s) is longer 
than L/2; likewise, the arclength s_(£/2) is shorter than 
L/2. In the region < s < L/2, the material points of 
r + have slid backwards relative to the material points of 
r_. Integrating ds+ — ds_, we find that the distance 
filament r_ slides past filament r + at s is 

A(s) = - [ aK(s)ds = -a[9(s)-9(0)}, (10) 
Jo 

where 9 is the angle between ei and e x . Note that 
Eq. (|10p is exact, although from here on we assume 
an <C 1. 

The composite filament is actuated by molecular mo- 
tors that slide r± back and forth past each other. Let 
/ m denote the force per unit length in the ei direction 
that the upper filament r + exerts on the lower filament 
r through the motors. To first order in deflection, the 
total moment acting across the cross-section at s is there- 
fore given by 

M z =Ah"(s) + a[ f m (s)ds, (11) 

J s 

In this model, the bending modulus A is an effective 
modulus representing the stiffness of the entire flagellar 
structure. We note that there are additional forces that 
may be included in modeling eukaryotic flagella, such as 
those arising from linking proteins that resist the slid- 
ing of filaments past each other. In this work we leave 
those out for simplicity; we have found that for reason- 
able magnitudes of these forces they do not change our 
results qualitatively. Moment balance, Eq. @, applies 
as before, leading to a force on the cross section at s of 
F = (—Ah'" + af m )e y . or an internal force per unit length 
oif int = (-Ah"" + afl n )e y . 

The boundary conditions for the active filament involve 
extra terms from the sliding force. At the free end, zero 
force and moment give 

-Ah"'(L)+af m (L) = (12) 
Ah"{L) = 0. (13) 

At the attachment to the head we may consider a variety 
of situations. For a fixed head, we may consider clamped 
and hinged boundary conditions. For clamped boundary 
conditions, 

h(0) = (14) 
h'(0) = 0. (15) 



For hinged boundary conditions where the total moment 
vanishes, 

Ah"(0) + a [ f m (s)ds = (16) 
Jo 

h{0) = 0. (17) 



B. Linearly and nonlinearly viscoelastic fluids 

In the limit of zero Reynolds number, the equation 
of motion of an incompressible medium is given by force 
balance, — Vp+V-T = 0, where pis the pressure enforcing 
incompressibility, V • v = 0; v is the velocity; and r is 
the deviatoric stress tensor. 

The material properties of the medium are given by 
the constitutive relation, which specifies the stress in the 
medium as a function of the strain and strain history of 
the material. In a Newtonian fluid with viscous forces, 
the stress tensor is proportional to the rate of strain. In 
contrast, in a viscoelastic fluid, there is an additional 
elastic component to the response of the medium. This 
effect is incorporated by fading-memory models, in which 
the stress relaxes over time to the viscous stress. We will 
focus on a particular fading-memory model, the Oldroyd- 
B fluid, which has the constitutive relation 

v v 
t + A t= 7/7 + Ar/ S 7 . (18) 

v 

In this equation, A is the single relaxation time, T= 

dtT + v- Vt — (Vv) T -t — t • Vv is the upper-convected 

time derivative of t, v is the velocity, rj is the polymer 

T v 
viscosity, 7 = Vv + (Vv) is the strain rate, and 7 is 

the upper-convected time derivative of the strain rate. 
In this paper we will use a dot to denote time differen- 
tiation. The total viscosity rj = rj p + r] s is the sum of 
polymer and solvent viscosities. The Oldroyd-B fluid is 
one of the simpler nonlinear constitutive relations for a 
viscoelastic fluid. The upper-convected time derivative 
is the source of nonlinearities If the solvent viscos- 
ity is ignored (r] 8 = ), the Oldroyd-B fluid reduces to 
the Upper Convected Maxwell fluid, and if there is no 
polymer (77 = rj s ) the Oldroyd-B fluid reduces to a New- 
tonian fluid. The Oldroyd B fluid is known to describe 
elastic and first normal stress effects as long as elonga- 
tional flows are not large; however, it does not capture 
shear-thinning or yield-stress behaviors observed in some 
non-Newtonian fluids. 

Although we have introduced the nonlinearly viscoelas- 
tic fluid, in what follows we will only work to linear order 
in the displacements of filaments. In that case, our re- 
sults will be the same as if we consider the same model 
with ordinary time derivatives, which corresponds to a 
linear Maxwell-Kelvin fluid. 
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III. FILAMENTS IN NEWTONIAN FLUIDS 

The motion and swimming properties of filaments in 
Newtonian fluids have been actively studied by many re- 
searchers @, B Eii El ■ Here we briefly summarize results 
in the Newtonian case which will be useful for comparison 
to filaments in viscoelastic fluids. 

A common approach to calculating the viscous forces 
acting on a flagellum is resistive force theory, a local drag 
theory in which the forces per unit length acting on a thin 
filament are proportional to the fluid velocity relative to 
the filament [20] ■ Although resistive force theory does not 
completely capture the effects of hydrodynamic interac- 
tions, macroscopic-scale experiments have shown that it 
is surprisingly accurate for single filaments [H, H3, I22I ]. 
For a purely viscous liquid, the force per unit length act- 
ing on a filament moving in an otherwise quiescent fluid 
is given by 

fvi s = -C||(r)||-C±(r)x, (19) 

the subscripts || and _L denote the directions along and 
perpendicular to the filament at r(s, t), respectively. The 
friction coefficients are proportional to the viscosity 77, 
with Cl « 2C||. 

For simplicity, assume that there is only one frequency 
of motion for the filament, so the displacement of the 
filament can be written h(s,t) = Re{/i(s, t) exp(— iu>t)}. 
In this case, using the resistive force theory, the y- 
component of the equation of motion for the passive and 
active filament are, respectively, 

-koCxh = -Ah"" (20) 
-iuC±h = -Ah"" + af! nt . (21) 

The solutions of these can be obtained in a manner sim- 
ilar to that described in Sec. |V]to give h(s, t). 

From h(s, i), we may calculate the power dissipated by 
the filament and the swimming velocity of the filament. 
At each material element of the filament, the power dissi- 
pated is the inner product of the hydrodynamic force and 
the velocity. The total time-averaged power dissipated is 
therefore 

(P) = ^f£ds\h(s)\ 2 (22) 

The swimming velocity can be determined by the con- 
straint under Stokes flow that the total force on the swim- 
mer is zero. To lowest order, the component of force in 
the e y direction is zero as a result of the equation of mo- 
tion for the filament shape. The beating motion of the 
filament produces a e x component of the time-averaged 
force that is second order in the deflection, which is bal- 
anced by a drag force from overall translation, i.e. swim- 
ming with velocity Z7rft i n the e x direction: 

C/rft = ^^vT^ Im [ ds h'h* (23) 
2iC|| Jo 



In this paper we ignore the drag from the head, which 
would contribute an additional factor in the expression 
for ZTrft- Note that for clamped boundary conditions, 
the torque exerted on the filament at s = has a reaction 
torque that causes the head and ultimately the whole 
sperm to rotate ■ Here we ignore these effects for two 
reasons: 1) Experiments looking at the shape of beating 
filaments can be performed with an immobilized head. 2) 
For swimming, calculations ignoring the head are simpler 
and can be easily compared to previous work, e.g. [Ill ] . 
and also are more easily extended and compared to the 
calculations we will perform in the viscoelastic fluid. 

It is important to note that it is only necessary to 
solve for the beating motion to first order in deflection 
amplitude to obtain the correct leading order results for 
the power and the swimming speed, even though these 
quantities are of second order in deflection. 



IV. FILAMENTS IN VISCOELASTIC FLUIDS 

Having seen how to analyze filaments in Newtonian 
fluids, we turn to viscoelastic fluids. In a viscoelastic 
fluid, the hydrodynamic force felt by the filament is no 
longer given by Eq. (TT9"]) . Due to the change in forces, 
the equations of motion for the filament are changed, so 
that the beating patterns change. Since both the beat- 
ing patterns and forces affect the power dissipated and 
swimming velocity, we expect both of those to change 
as well. To calculate the power dissipated by a filament 
beating in a viscoelastic medium to leading order in the 
deflection, we only need the velocity and the force to 
first order. Thus, it is sufficient to find the shape to 
first order. From Eq. (|23[) we expect that a change in 
the first-order shape h will change the swimming veloc- 
ity; however, there are additional second order correc- 
tions that are present even for a waveform of prescribed 
shape [HI, [l7|. Consistency requires including both of 
these second order corrections [la ]. In this paper we 
concentrate on the effects of viscoelasticity on beating 
patterns and power dissipated. 

The calculation of the force per unit length acting on 
a slender filament undergoing large deflections in a vis- 
coelastic fluid is a daunting challenge. Since we consider 
small deflections, we may consider the simpler problem of 
the forces per unit length acting on a slightly perturbed 
infinitely long cylinder with a prescribed lateral or longi- 
tudinal traveling wave. Elsewhere we solve for the flow 
induced by these motions using the Oldroyd-B fluid and 
calculate the forces acting on the deformed cylinder [23| . 
We find that the force per unit length foidroyd obeys 

fbhoyd + Al foidroyd = fvis + — Afvis. (24) 

In the limit of zero solvent viscosity, r] s = 0, Eq. ([24]) is 
the same as the resistive force theory proposed by Ful- 
ford, Katz, and Powell for the linear Maxwell model [lj|- 
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Eq. [24] is valid for slender filaments and to first order in 
small deflections. 



V. RESULTS 
A. Non-dimensionalization 



A. Application to passive and active filaments 

To apply the above results to the passive and active 
filament models, we use the hydrodynamic force foidroyd 
in the equations of motion. For the passive filament with 
a single frequency of motion, we find 



. l-iDe 2 ,, 

ilu Cih = —An 

1 -iDe S 



(25) 



where the Deborah numbers De = lo\ and De2 = Derj s /r). 
For the active filament we must include the motor forcing: 



f - i De 2 

I ■ 

1 -iDe 



Uh = -Ah"" + af! n 



(26) 



The beating pattern is determined by solving these equa- 
tions of motion, and once the beating pattern is known, 
the power dissipated and swimming velocity can be cal- 
culated. 

The power dissipated is calculated in the same way as 
for a Newtonian fluid, except that now the force per unit 
length is specified by Eq. (f2~4"|) , leading to 



(P) 



C±w 2 1 + DeDe 2 
~~ 2 1 + De 2 



ds \h(s)\ 



(27) 



As noted above,once the beating pattern affects the 
swimming velocity. However, there is an additional cor- 
rection to the swimming velocity in the nonlinearly vis- 
coelastic Oldroyd-B fluid. Instead of Eq. [23l the correct 
expression for the swimming speed is [15l . l23j 



U 



1 + DeDe 



1 + De \j 
1 + DeDe 2 
1 + De r 



ds h'(s,t)h(s,t)J (28) 
j U RFT , (29) 



where the second line follows from the ratio Ci/£|| = 2. 
This result applies to an infinite cylinder moving with 
beating pattern h(s,t) = Re{/i(s) exp(— iuit)}. To apply 
this result to the finite filament, we imagine periodically 
replicating the calculated beating pattern of the filament 
so that it becomes infinitely long. Note that due to the 
periodicity, the integration only needs to be performed 
from to L. This result ignores end effects from the flow 
around the end of the finite filament. We have mentioned 
the effects of beating patterns on swimming velocity here 
for completeness. In the rest of the paper we focus on 
the changes to beating patterns themselves and power 
dissipation rather than swimming velocity. 



We can non-dimensionalize the filament equations of 
motion and boundary conditions by measuring lengths in 
terms of L, f m in terms of A/aL 2 , and time in terms of 
uj^ 1 . Velocities are measured in terms of Lui and angular 
velocities in terms of u. For notational simplicity, after 
scaling we use the same symbols for the new quantities. 
The non-dimensional equation of motion for the passive 
filament with a single frequency is 



1 — iDe 
with boundary conditions 

-h"'(l) = 
h"(l) = 
h(0) = 



(30) 



(31) 
(32) 
(33) 



and either of 

h'(Q) 
h"'(0) 



ecos(i) (prescribed angle) (34) 
ecos(i) (prescribed force) . (35) 



The dimensionless parameter Sp = L (lo^/A) 1 ^ 4 , known 
as the Sperm number, involves the ratio of the bending 
relaxation time of the filament to the period of the trav- 
eling wave. In what follows the forcing e is chosen small 
enough so that we remain safely in the linear regime. 

The non-dimensional equation of motion for the active 
filament with a fixed head and single frequency becomes 



4 l-iDe 2 ~ -„„ 

iSp ——h + h 

1 — iDe 



f = 



(36) 



where f m = Re{/ m (s) exp(— it)}. The boundary condi- 
tions for Eq. (|36p are 



h"'(l) + UD 

h"(l) 

h(0) 



and either of 



h"(Q) 



h'(0) = (clamped) 
fm(s)ds = (hinged). 



(37) 
(38) 
(39) 

(40) 
(41) 



B. Passive filament 



For the passive filament in sinusoidal steady state, we 
examine the solutions to Eq. <PU|) which are 



4 
3=1 



(42) 
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J3e 



20 



100 



6= 0.2 



6= 0.1 




e = 1.0 
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FIG. 3: (Color online.) Shapes of beating patterns for passive 
filaments driven by prescribed angle at the left end (s — 0). A 
half-cycle (red, orange, green, light blue, blue) of the pattern 
is shown for viscous (De = 0) and viscoelastic (De = 100) 
cases, and Sperm number 2, 7, and 20. The ratio of solvent 
viscosity to the total viscosity is r) s /r) — De2/De = 10 -4 . At 
the top of each plot, the (dimensionless) angular magnitude 
e required to produce motion with amplitude 0.1L is shown. 



where the ?%■ solve 



iSp 



1 -iDc 2 
1 -iDe 



0. 



(43) 



The four coefficients Cj are determined by the boundary 
conditions, and the full time-dependent solution is then 
given by 



h{s,t) = Re{e~ i */i h (s)} 



(44) 



From the form of hh in Eq. (|42|) , we can already make a 
general comment about the expected form of beating pat- 
terns. The exponentials in Eq. (j4"2"|) will lead to traveling 
or standing waves. Standing waves can only be produced 
if there are pairs of Tj that are complex conjugates. This 
happens only in the elastic limit (De — > oo) for which the 
Tj are proportional to ±(1 ± i). Therefore in the viscous 
case the beating pattern will always have the form of a 
traveling wave, but as De increases the beating pattern 
takes on more characteristics of a standing wave. 

Representative plots of these solutions for prescribed 
attachment angles and forces are presented in Fig. [3] 
Solutions for the viscous (De = 0) and viscoelastic 
(De = 100) limits are plotted side by side. As ex- 
pected, changing the viscoelastic character of the fluid 
by varying the Deborah number induces changes in the 
shapes of the beating patterns. An important feature of 



these solutions is the presence of a bending length scale 
£ £/L = \rj\- 1 = Sp|(l - iDe 2 )/(l - iDe)!" 1 / 4 . For 
De = 0, the purely viscous case, this corresponds to the 
length scale set by the Sperm number; a stiff filament 
with small Sperm number has large £ and acts as a rigid 
rod, while a flexible filament with large Sperm number 
has small £ and can bend in response to viscous forces. 
Increasing the Deborah number acts to increase £. In the 
figure, the amplitude of the prescribed driving angle is se- 
lected so that the maximum displacement of the filament 
is L/10, and then printed above each figure. The beating 
pattern is linear in the amplitude of the driving angle. 
Due to larger bending lengths, £, in viscoelastic fluids, 
smaller driving angles are typically required to produce 
the same amplitude motion. 

The beating patterns in Fig. [3] also demonstrate that 
viscoelastic effects alter the nature of the beating pat- 
terns as expected from the form of Eq. I|42p: for De = 0, 
we obtain traveling wave beating patterns, while for large 
De, in the elastic limit, we obtain standing wave beating 
patterns. 

From Eq. (|27[) we obtain the time-averaged power dis- 
sipated 



(Pi 



Ml 



Cj 



Re 



dsuj \h\ 



1 - i Deo 



iDe 



1 + DeDe 2 
2(1 + De 2 ) 



(45) 
(46) 



Figure 2] shows the power dissipation versus Deborah 
number (dashed curve). We also display the power dis- 
sipation, Eq. (|46p . evaluated with the purely viscous 
(De = 0) shape h (solid curve). Note that the difference 
in shape h between the viscous and viscoelastic cases has 
little effect; the power is mainly governed directly by the 
dependence of the force per unit length on De. 



C. Active Filaments 

1. Fixed Head 

Examining Eq. (|36[) . we find that it only differs from 
the passive case at the boundary conditions and by the 
presence of an inhomogeneous term. Therefore we may 
use the homogeneous solution from the previous subsec- 
tion, 



hh(s) 



4 

£ 

3=1 



C l exp(r,-s) 



(47) 



and add to it a particular solution. As in the passive 
case, beating patterns with traveling waves are expected 
in the viscous limit, and increasing the Deborah number 
leads to beating patterns with characteristics of standing 
waves. 
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FIG. 4: Time- averaged power dissipated by passive filaments 
driven by prescribed angle, as a function of A = lu\ (dashed 
curve). Power is normalized to the power dissipated in the 
viscous case De = 0. The black curve corresponds to the 
power dissipated if the beating pattern does not change from 
the viscous case, but the filament is placed in a viscoelastic 
medium. The driving amplitude is chosen so that the maxi- 
mum amplitude of the beating pattern is L/10 in the viscous 
case. The ratio of solvent viscosity to the total viscosity is 
Vs/v = De 2 /De = 1(T 4 . 

To find the particular solution, let us assume that the 
internal sliding force takes the form of a sliding wave 

f m (s) = Re {/ m exp (iks - it)} (48) 

In that case, we can write the particular solution 

hp = ifc/ m e ife 7 [-iS P 4 (l -iDe 2 )/(l -iDe) + k 4 ] . 

(49) 

Again the four boundary conditions set the four unknown 
coefficients and the full time-dependent solution is 

h(s,t) =Re{e- it [h h (s) + h p ]} (50) 

The beating shapes of the active flagellum and a fixed 
head are shown in Fig. [5j for De = (viscous case) and 
De = 100 (viscoelastic case), and for various spatial de- 
pendencies of the internal sliding force. The forces re- 
quired to produce beating with amplitude 0.1L are larger 
in the viscous cases than the viscoelastic cases, again 
since the bending length £ is longer in the viscoelas- 
tic case than the viscous case. Just as in the case of 
the passive filament, the effect of the length scale £ can 
be seen by comparing the viscous and viscoelastic cases. 
However, the wavelength 2ir/k of the sliding force intro- 
duces another length scale. The combined presence of 
two different length scales can be seen, for example, in 
the shape of the hinged filament with De = 0, and (di- 
mensional) k — 8n/L. In this example, £ is longer (and is 
the same as that seen for all the Dc = shapes of Fig.[SJ), 
while the shorter scale 2ir/k leads to ripples on top of the 
longer deformation. Comparing the viscous to viscoelas- 
tic case, one can also see a decrease in the visible effects 
of the shorter lengthscale arising from the internal forces 
as £ increases and the shape is dominated by drag forces. 



For example, consider Fig. [6l which shows the effect of 
varying the Sperm number for an internal sliding force 
with (dimensional) k — 8tt/L. As the Sperm number 
increases, £ decreases and smaller lengthscale motions 
become more visible. One implication of the presence of 
these two length scales is that the observed "wavelength" 
of a beating pattern does not necessarily give direct in- 
formation about either lengthscale. 

The time-averaged power dissipated by the clamped 
and hinged internally driven filaments can be calculated 
as in the previous section. The power versus De, normal- 
ized by viscous power with A = 0, is shown in Fig. [7k 
and [7Jd. Also shown for comparison is the power dis- 
sipated assuming that the beating shape of the viscous 
case is prescribed for the viscoelastic cases. For both 
prescribed forces and velocities, the growing imaginary 
component of the viscosity shifts the drag force to be out 
of phase with the velocity, decreasing the power dissi- 
pated. However, for large De the power dissipated for 
prescribed forces is always greater than that dissipated 
for prescribed velocities, because as De increases, the real 
part of the drag decreases in amplitude. Thus, prescribed 
forces lead to motion with larger amplitude and velocity, 
and therefore more power dissipated relative to the case 
with prescribed velocities. 



2. Comparison to experimental beating patterns 

The equation of motion for the active filament (Eq. [36]) 
allows us to understand some of the qualitative changes in 
the shapes of beating patterns of sperm flagella in differ- 
ent media. There have been a number of studies observ- 
ing the different beating patterns of sperm in media with 
increased viscosity and viscoelasticity [3, Of these, 
we focus on the paper of Ishijima et al., because it com- 
bines a systematic study of beating patterns in different 
media with measurements of viscosity using a falling-ball 
rheometer, and because its use of a pipette to hold the 
heads of observed sperm still makes its results amenable 
for modeling. In Fig. [8k, the beating patterns of human 
sperm in Hanks' medium with viscosity of 1 cp, 35 cp, 
4000 cp, and cervical mucus with viscosity 4360 cp are 
shown. The 35 cP and 4000 cP viscosity media are made 
by adding polyvinylpyrrolidone and methylcellulose, re- 
spectively, to Hanks' medium (a Newtonian solution of 
salts and glucose in water). As viscosity increases, mo- 
tion is concentrated at the distal end, and motion in the 
proximal and middle portions of the flagella is damped 
out. Similar behavior was observed in Refs. UHI. 

At very high viscosity, asymptotic analysis of the equa- 
tion of motion Eq. ([3"6"| demonstrates that most of the 
beating motion occurs at the end of the flagellum. For 
concreteness, because the beating patterns in Fig. [5k 
were obtained while the sperm head was held in place by 
a pipette, we use clamped boundary conditions to model 
the motion. Furthermore, since wavelengths of 26-28 /im 
were observed in the low viscosity media for flagella of 



8 



Del 

k/b 



371 



871 



f„ = 119.0 



Clamped 







s/L 



100 





£ m = 4.3 




. 



f m = 14.0 



f m = 16.4 



f m = 123.2 



Hinged 







100 
























fm = 


2 . 3 


1 



f.= 21.5 





FIG. 5: (Color online.) Shapes of beating patterns for filaments with internal sliding forces and fixed head position, with Sp = 7. 
A half-cycle (red, orange, green, light blue, blue) of the pattern is shown for viscous (De = 0) and viscoelastic (De = 100) cases, 
and for internal sliding forces given by Eq. (|48|l with k/L varying from (uniform force) to 8n, as indicated. At the top of 
each figure, the (dimensionless) magnitude / m required to produce motion with amplitude 0.1L is shown. The ratio of solvent 
viscosity to the total viscosity is rj a /r] — De2/De = 10 -4 . 



overall length « 40/im, we assume a sliding force with 
(nondimensional) k = 3tt and the form of a traveling 
wave 

f m (s)=Re{f m e iks - it } (51) 

Then the inhomogeneous part of Eq. (|3"6"|) has the partic- 
ular solution 



hp = ifc/ m e Ifc 7 [-iS P 4 (l - iDe 2 )/(l -iDe) + fc 4 ] 
while the homogeneous part hh(s) obeys 



4 1 - i De 2 
l Sp — , ^ /ih + h 



mi 
h 



1 - iDe 
with boundary conditions 

= / m (i)-<(i) 
^(i) = -^'(l) 

h h (0) = -h p (Q) 
h' h (0) = -h' p (0). 



(52) 
(53) 



(54) 
(55) 
(56) 
(57) 



We analyze the homogeneous solution by asymptotically 
matching solutions appropriate for the middle of the flag- 
ellum with solutions appropriate for "boundary layers" at 
the ends of the fiagcllum. At high viscosity, the sperm 
number is large, and therefore the term involving the 
fourth power of the sperm number in Eq. (|53[) dominates. 
Throwing out the other terms, we find that the homoge- 
neous solution vanishes. This balance is valid for the 
middle of the fiagellum, but in a boundary layer near 
the ends of the fiagellum, the homogeneous solution may 
vary rapidly, making the term with derivatives impor- 
tant. The size of this boundary layer can be expected to 
be of order £ = L\(l - iDe)(l - iDe 2 )| 1/4 /Sp. Measur- 
ing lengths by £, i.e. x = s/£ for the proximal region, 
x = (L— for the distal region, we obtain the equation 
of motion 



c = 



-e ie h h + diK 
1-iDea 



1 -iDe 



1 - iDe 



1 - i De 2 



(58) 
(59) 
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FIG. 6: (Color online.) Shapes of beating patterns for hinged 
filaments with internal sliding forces, with varying Sp. A 
half-cycle (red, orange, green, light blue, blue) of the pattern 
is shown for viscous (De = 0) and viscoelastic (De = 100) 
cases, and for internal sliding forces given by Eq. [48] with 
k/L — 8n. At the top of each figure, the (dimensionless) 
magnitude f m required to produce motion with amplitude 
0.1L is shown. The ratio of solvent viscosity to the total 
viscosity is i] s /r] = De2/De = 10~ 4 . 



At the proximal end, the boundary conditions are 



K{x = 0) 
d x h h (x = 0) 
hh(x — > oo) 



-Mo) 

-d x h p (x 



0) 







(60) 
(61) 
(62) 



While at the distal end the boundary conditions are 

dlK{x = 0) = £ 3 / m (l) - d 3 x h p (x = 0) (63) 



5gfc h (s = 0) = -Sgftp(as = 0) 
hh{x — > oo) = 



(64) 
(65) 



In the boundary conditions, the last constraint comes 
from the requirement that the solution in the boundary 
layers asymptotically matches the solution in the middle 
of the flagellum. The solution to these equations is a 
linear combination of four exponential functions. 



K(x) = ^2 c„exp(p„x) 

n=0,l,2,3 

p n = exp(i#/4 + irm/2) 



(66) 
(67) 



The p n have real and imaginary parts of similar order, 
so that the solutions decay or increase and have spatial 
oscillations on a similar lengthscale, which is approxi- 
mately £ in dimensional units. Due to the asymptotic 



FIG. 7: Time-averaged power dissipated by filaments with ac- 
tive internal sliding forces, as a function of De = lo\. Power is 
normalized to the power dissipated in the viscous case De = 0, 
with force amplitude such that maximum filament displace- 
ment in the viscous case De = is 0.1L. Different curves 
correspond to sliding forces given by Eq. [48] with different 
values of k: dashes, k = 0; long dashes, k — tt/L; long dash- 
short dash, k = 3n/L. The solid curve is the power dissi- 
pated by a filament with a beating shape corresponding to 
the viscous case (De = 0). The ratio of solvent viscosity to 
the total viscosity is r/ s /ri = De2/De = 10 -4 . a) Fixed head 
with clamped attachment point, b) Fixed head with hinged 
attachment point. 



matching conditions, only the exponentials that decrease 
as x increase are allowed, corresponding to po and p3 if 
< 9 < 2tt. 

At the proximal end, due to the clamped boundary 
conditions the homogeneous solution must cancel the in- 
homogeneous solution at the attachment to the head. 
The homogeneous solution is specified by 



Co 



C3 



-M P (o) + ifcM°) 

P3 - P0 

p o h p (0) - ifcfop(O) 

P3 - P0 



(68) 
(69) 



which has an amplitude of order / m (£/£) 4 , the same 
as the inhomogeneous solution. At the distal end, the 
largest term in the boundary conditions comes from bal- 
ancing the sliding force. The homogeneous solution is 
specified by 



Co 



-p 3 (ifc) 2 fe P (i) - leup - (ifc) 3 yi)] 

Pl(P3 - Po) 

Po (ifc) 2 fe p (i) + K 3 / m (i)-(ifc) 3 fep(i)] 

Pl(P3 - Po) 



(70) 
(71) 
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FIG. 8: (Color online.) a) Beating patterns of human sperm observed by Ishijima et al. [j, in (from left to right) lcP Hanks' 
solution, 35 cP Hanks' solution, 4000 cP Hanks' solution, 4360 cP cervical mucus, b) Shapes of beating patterns for filaments 
with internal sliding forces, fixed head position, and clamped boundary conditions. A half-cycle (red, orange, green, light blue, 
blue) of the pattern is shown for internal sliding forces given by Eq. [51] with k/L — 3n From left to right, the sperm numbers 
are 5.6, 13.6, 38.7, and 45, corresponding to viscosities of 1 cP, 35 cP, 4000 cP, and 4360 cP. In the first two plots, the medium 
is purely viscous, while in the third plot, De = 0.5 and De2 = 0.5/4000. In the fourth plot, De = 50 and De2 = 50/4360. At 
the top of each figure, the (dimensionless) magnitude f m required to produce motion with amplitude 0.1L, and the lengthscale 
£ is shown. 



which has amplitude of order / m (£/L) 3 . For high viscosi- 
ties and low to moderate Deborah number this amplitude 
can be much larger than the amplitude of the proximal 
and inhomogeneous solutions. In the viscous case, with 
De = 0, and high viscosity, the ratio between the ampli- 
tude of the homogeneous and inhomogeneous solution in 
the distal portion of the flagellum is of order fc/Sp. 

To summarize, in the middle of the flagellum the ho- 
mogeneous solution vanishes and the beating shape is 
dominated by the inhomogeneous solution, which has a 
magnitude which decreases as viscosity increases, and is 
quite small for high viscosity. At the proximal end of 
the flagellum, the homogeneous solution is of the same 
order as the inhomogeneous solution due to the clamped 
boundary conditions. On the other hand, the motion at 
the distal end is dominated by the homogeneous solution, 
has larger amplitude than the middle portion, decays ex- 
ponentially with a lengthscale £, and can have oscillations 
with wavelengths up to a few times smaller than £. Thus, 
in high viscosity solutions the motion of flagellar beating 
patterns are constrained to the distal tip. Although we 
assume only a single mode of the sliding force, the asymp- 
totic analysis is generally valid, and should also apply to 
more realistic sliding forces. 

In Fig. [SJd we plot flagellar beating patterns for a slid- 
ing force of k — 3ir/L. We constrain the motion to 
lie in a two-dimensional plane for simplicity, although 
the observed beating patterns are three-dimensional. 
These beating patterns are obtained by solving the linear 
Eq. exactly, i.e., not using the asymptotic analysis. 
In this modeling we assume the pure Hanks' solution and 
Hanks' solution with polyvinylpyrrolidone to be Newto- 
nian. For the Hanks' solution with methylcellulose we as- 
sume a Deborah number of 0.5, corresponding to a beat- 
ing frequency of 7 Hz [J] and a relaxation time constant 



of about 10~ 2 s 24]. For the cervical mucus we assume 
a Deborah number of 50, corresponding to a beating fre- 
quency of 12 Hz [4| and a relaxation time constant of a 
little less than a second [25[ . To determine De2 we use the 
ratio rjs/i] =1 cP/4360cP. The Sperm numbers are deter- 
mined from L = 40 /mi, C_L = 2r], A = 4x 10~ 22 Nm 2 [ill ], 
viscosities from Fig. [5J and frequencies of 12 Hz for the 
lcP and 35 cP media, 7 Hz for the 4000 cP media, and 
12 Hz for cervical mucus [4|. 

In the plots, the amplitude of the beating pattern is 
proportional to the driving force / m . The amplitude in 
the middle sections is therefore nearly completely sup- 
pressed as viscosity increases as compared to the ampli- 
tude in the 1 cP case. At the same time, the motion near 
the end of the flagellum is relatively unsuppresed. In ad- 
dition the observed beating pattern in cervical mucus has 
motion extending farther away from the distal end than 
the beating pattern in 4000 cP Hanks' solution, in both 
the experimental (Fig. [HJl) and modeled (Fig. [5Jd) beating 
patterns. This agrees with the fact that the lengthscale 
£ decreases as viscoelastic effects become important. 



VI. DISCUSSION 

In a viscoelastic medium the forces exerted on a flexi- 
ble swimmer are different from those exerted by a Newto- 
nian medium. We have calculated the viscoelastic forces 
for a medium with fading memory using resistive force 
theory. The effects of hydrodynamic forces on flexible 
swimmers can be seen in the passive filament model we 
have described. We have also used a simple model of an 
active filament to shed some light on the motion of sperm 
flagella in different media. 

In agreement with experimental observations, our slid- 
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ing filament model shows that the beating is confined 
to the distal tip for the high viscosity Hanks' solution. 
While the experimental beating patterns are qualitatively 
explained by the high- viscosity asymptotic analysis, and 
the solutions in Fig. [8)3 show the same qualitative trends, 
our model is not quantitatively accurate. For example, in 
the cervical mucus, our model shows much less confine- 
ment of motion to the distal tip than the observed beat- 
ing patterns. In particular, we note that according to our 
modeling, only the third panel of Fig.[HjD is in the strongly 
asymptotic regime, in contrast to the beating patterns of 
Ishijima et al, in which the motion in cervical mucus also 
seems to be in the asymptotic regime. The discrepancy 
between our models and the observed beating patterns 
may be due to an overestimation of the Deborah number, 
our approximation of small amplitudes (the actual beat- 
ing patterns show rather large curvatures), or our use of 
a simple form of the sliding force. Ishijima et. al. observe 
varying beating frequencies in different media, indicating 
that the sliding force itself is changing in response to the 
different drag forces exerted in different media. Changes 
in the active force due to changes in the media have been 
modeled by various workers [HI, [2(| [27j ■ We have spec- 
ified a fixed force; however, our analysis might be useful 
in extending studies of the sliding mechanism such as in 
Ref. [IH to high viscosity and viscoelastic situations. Fi- 
nally, our simple model does not take into account the 



possibility of non-Newtonian behavior which extends be- 
yond viscoelasticity in Hanks' medium with added poly- 
mers. We note that in the case of methylcellulose solu- 
tions, the swimming of bacteria with helical flagella has 
been analyzed by introducing different viscosities for par- 
allel and perpendicular drag coefficients [ll]- This type 
of treatment would not affect our results, since at first 
order the shape is determined solely by the perpendic- 
ular drag component in our model. However, the work 
of reference [28| points out that care should be taken to 
measure appropriate viscosities for flagellar movements, 
especially as the viscosity of methycellulose solutions is 
known to be shear-rate-dependent. 

Our work is a first step to understanding how swim- 
ming is modified in a complex medium. There is scope to 
expand these studies, probably numerically, in addressing 
questions such as the role of end effects and associated 
elongational flows, the role of drag from a head and the 
role of large amplitude motion in viscolastic fluids. In 
addition, many biological media, including mucus in the 
reproductive tract, are in gels, which may be expected to 
have different effects on the swimming than fluids. 
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and DMS-0615919 (TRP), and NIH grant No. R01 
GM072004 (CWW). TRP and HCF thank the Aspen 
Center for Physics, where some of this work was done. 
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